#include<cstdlib>
#include<cmath>
#include<exception>
#include<stdexcept>
#include<string>
#include<iostream>
#include<fstream>
#include<vector>
#include<numeric>

#include<boost/format.hpp>

#include<formic/lapack/interface.h>

namespace formic {

  namespace test {

    template <class S> inline void equal(const S & a, const S & b, const std::string & msg) {
      if ( double(std::abs(a-b)) < 1.0e-9 )
        return;
      std::string error_str("formic::test::equal failed for ");
      error_str.append(msg);
      throw std::runtime_error(error_str);
    }

    template <class S> inline void equal(const std::vector<S> & a, const std::vector<S> & b, const std::string & msg) {
      std::string error_str("formic::test::equal failed for ");
      error_str.append(msg);
      if ( a.size() != b.size() )
        throw std::runtime_error(error_str);
      for (size_t i = 0; i < a.size(); i++)
        if ( double(std::abs(a.at(i)-b.at(i))) >= 1.0e-9 )
          throw std::runtime_error(error_str);
    }

  }

}

void formic_lapack_tests() {

  // dscal
  {
    const int n = 4;
    std::vector<double> v0(n, 1.0);
    formic::dscal(n, 0.5, &v0.at(0), 1);
    formic::test::equal(2.0, std::accumulate(v0.begin(), v0.end(), 0.0), "dscal");
    formic::xscal(n, 0.5, &v0.at(0), 1);
    formic::test::equal(1.0, std::accumulate(v0.begin(), v0.end(), 0.0), "xscal, x=d");
  }

  // zscal
  {
    const int n = 4;
    const std::complex<double> zero(0.0, 0.0);
    const std::complex<double> answer(4.0, 7.0);
    std::vector<std::complex<double> > v0(n, std::complex<double>(2.0, 1.0));
    formic::zscal(n, std::complex<double>(0.75, 0.5), &v0.at(0), 1);
    formic::test::equal(answer, std::accumulate(v0.begin(), v0.end(), zero), "zscal");
    std::vector<std::complex<double> > v1(n, std::complex<double>(2.0, 1.0));
    formic::xscal(n, std::complex<double>(0.75, 0.5), &v1.at(0), 1);
    formic::test::equal(answer, std::accumulate(v1.begin(), v1.end(), zero), "xscal, x=z");
  }

  // dcopy
  {
    const int n = 4;
    std::vector<double> v0(n, 1.5);
    std::vector<double> v1(n, 0.0);
    std::vector<double> v2(n, 0.0);
    formic::dcopy(n, &v0.at(0), 1, &v1.at(0), 1);
    formic::test::equal(6.0, std::accumulate(v1.begin(), v1.end(), 0.0), "dcopy");
    formic::xcopy(n, &v0.at(0), 1, &v2.at(0), 1);
    formic::test::equal(6.0, std::accumulate(v2.begin(), v2.end(), 0.0), "xcopy, x=d");
  }

  // zcopy
  {
    const int n = 4;
    const std::complex<double> zero(0.0, 0.0);
    const std::complex<double> answer(6.0, 1.2);
    std::vector<std::complex<double> > v0(n, std::complex<double>(1.5, 0.3));
    std::vector<std::complex<double> > v1(n, std::complex<double>(0.0, 0.0));
    std::vector<std::complex<double> > v2(n, std::complex<double>(0.0, 0.0));
    formic::zcopy(n, &v0.at(0), 1, &v1.at(0), 1);
    formic::test::equal(answer, std::accumulate(v1.begin(), v1.end(), zero), "zcopy");
    formic::zcopy(n, &v0.at(0), 1, &v2.at(0), 1);
    formic::test::equal(answer, std::accumulate(v2.begin(), v2.end(), zero), "xcopy, x=z");
  }

  // daxpy
  {
    const int n = 4;
    std::vector<double> v0(n, 1.2);
    std::vector<double> v1(n, 0.1);
    std::vector<double> v2(n, 0.1);
    formic::daxpy(n, 2.0, &v0.at(0), 1, &v1.at(0), 1);
    formic::test::equal(10.0, std::accumulate(v1.begin(), v1.end(), 0.0), "daxpy");
    formic::xaxpy(n, 2.0, &v0.at(0), 1, &v2.at(0), 1);
    formic::test::equal(10.0, std::accumulate(v2.begin(), v2.end(), 0.0), "xaxpy, x=d");
  }

  // zaxpy
  {
    const int n = 4;
    const std::complex<double> zero(0.0, 0.0);
    const std::complex<double> answer(-0.4, 9.6);
    const std::complex<double> factor(0.5, 2.0);
    std::vector<std::complex<double> > v0(n, std::complex<double>(1.2,  0.4));
    std::vector<std::complex<double> > v1(n, std::complex<double>(0.1, -0.2));
    std::vector<std::complex<double> > v2(n, std::complex<double>(0.1, -0.2));
    formic::zaxpy(n, factor, &v0.at(0), 1, &v1.at(0), 1);
    formic::test::equal(answer, std::accumulate(v1.begin(), v1.end(), zero), "zaxpy");
    formic::xaxpy(n, factor, &v0.at(0), 1, &v2.at(0), 1);
    formic::test::equal(answer, std::accumulate(v2.begin(), v2.end(), zero), "xaxpy, x=z");
  }

  // ddot
  {
    const int n = 3;
    std::vector<double> v0(n);  v0.at(0) = 0.1;  v0.at(1) = 2.7;  v0.at(2) =  1.2;
    std::vector<double> v1(n);  v1.at(0) = 2.0;  v1.at(1) = 1.4;  v1.at(2) = -0.5;
    formic::test::equal(3.38, formic::ddot(n, &v0.at(0), 1, &v1.at(0), 1), "ddot");
    formic::test::equal(3.38, formic::xdot(n, &v0.at(0), 1, &v1.at(0), 1), "xdot, x=d");
  }

  // zdot
  {
    const int n = 3;
    const std::complex<double> z0( 1.2, -0.6);
    const std::complex<double> z1( 2.2,  1.6);
    const std::complex<double> z2(-0.3,  0.2);
    const std::complex<double> w0( 0.1, -0.1);
    const std::complex<double> w1( 1.4, -0.2);
    const std::complex<double> w2( 1.0,  0.5);
    std::vector<std::complex<double> > v0(n);  v0.at(0) = z0;  v0.at(1) = z1;  v0.at(2) = z2;
    std::vector<std::complex<double> > v1(n);  v1.at(0) = w0;  v1.at(1) = w1;  v1.at(2) = w2;
    const std::complex<double> answer(3.06, 1.67);
    formic::test::equal(answer, formic::zdot(n, &v0.at(0), 1, &v1.at(0), 1), "zdot");
    formic::test::equal(answer, formic::xdot(n, &v0.at(0), 1, &v1.at(0), 1), "xdot, x=z");
  }

  // ddotc
  {
    const int n = 3;
    std::vector<double> v0(n);  v0.at(0) = 0.1;  v0.at(1) = 2.7;  v0.at(2) =  1.2;
    std::vector<double> v1(n);  v1.at(0) = 2.0;  v1.at(1) = 1.4;  v1.at(2) = -0.5;
    formic::test::equal(3.38, formic::ddotc(n, &v0.at(0), 1, &v1.at(0), 1), "ddotc");
    formic::test::equal(3.38, formic::xdotc(n, &v0.at(0), 1, &v1.at(0), 1), "xdotc, x=d");
  }

  // zdotc
  {
    const int n = 3;
    const std::complex<double> z0( 1.2, -0.6);
    const std::complex<double> z1( 2.2,  1.6);
    const std::complex<double> z2(-0.3,  0.2);
    const std::complex<double> w0( 0.1, -0.1);
    const std::complex<double> w1( 1.4, -0.2);
    const std::complex<double> w2( 1.0,  0.5);
    std::vector<std::complex<double> > v0(n);  v0.at(0) = z0;  v0.at(1) = z1;  v0.at(2) = z2;
    std::vector<std::complex<double> > v1(n);  v1.at(0) = w0;  v1.at(1) = w1;  v1.at(2) = w2;
    const std::complex<double> answer(2.74, -3.09);
    formic::test::equal(answer, formic::zdotc(n, &v0.at(0), 1, &v1.at(0), 1), "zdotc");
    formic::test::equal(answer, formic::xdotc(n, &v0.at(0), 1, &v1.at(0), 1), "xdotc, x=z");
  }

  // create matrices using fortran indexing
  std::vector<double> md0(4);
  md0.at(0*2+0) = 3.0;
  md0.at(1*2+0) = 1.0;
  md0.at(0*2+1) = 2.0;
  md0.at(1*2+1) = 4.0;
  std::vector<double> md1(4);
  md1.at(0*2+0) = 1.0;
  md1.at(1*2+0) = 2.0;
  md1.at(0*2+1) = 2.0;
  md1.at(1*2+1) = 3.0;
  std::vector<double> md2(4);
  md2.at(0*2+0) = 7.0;
  md2.at(1*2+0) = 1.0;
  md2.at(0*2+1) = 1.0;
  md2.at(1*2+1) = 8.0;
  std::vector<std::complex<double> > mc0(4);
  mc0.at(0*2+0) = std::complex<double>( 3.0,  0.0);
  mc0.at(1*2+0) = std::complex<double>( 0.0,  1.0);
  mc0.at(0*2+1) = std::complex<double>( 0.0,  2.0);
  mc0.at(1*2+1) = std::complex<double>( 4.0,  0.0);
  std::vector<std::complex<double> > mc1(4);
  mc1.at(0*2+0) = std::complex<double>( 1.0,  0.0);
  mc1.at(1*2+0) = std::complex<double>( 0.0,  2.0);
  mc1.at(0*2+1) = std::complex<double>( 0.0, -2.0);
  mc1.at(1*2+1) = std::complex<double>( 3.0,  0.0);
  std::vector<std::complex<double> > mc2(4);
  mc2.at(0*2+0) = std::complex<double>( 7.0,  0.0);
  mc2.at(1*2+0) = std::complex<double>( 0.1, -0.2);
  mc2.at(0*2+1) = std::complex<double>( 0.1,  0.2);
  mc2.at(1*2+1) = std::complex<double>( 8.0,  0.0);

  // dgemm
  {
    const int n = 2;
    std::vector<double> mat(n*n, 1.0);
    formic::dgemm('N', 'N', n, n, n, 0.5, &md0.at(0), n, &md1.at(0), n, 0.7, &mat.at(0), n);
    formic::test::equal(3.2, mat.at(0*2+0), "dgemm");
    formic::test::equal(5.2, mat.at(1*2+0), "dgemm");
    formic::test::equal(5.7, mat.at(0*2+1), "dgemm");
    formic::test::equal(8.7, mat.at(1*2+1), "dgemm");
  }

  // xgemm, x=d
  {
    const int n = 2;
    std::vector<double> mat(n*n, 1.0);
    formic::xgemm('N', 'N', n, n, n, 0.5, &md0.at(0), n, &md1.at(0), n, 0.7, &mat.at(0), n);
    formic::test::equal(3.2, mat.at(0*2+0), "xgemm, x=d");
    formic::test::equal(5.2, mat.at(1*2+0), "xgemm, x=d");
    formic::test::equal(5.7, mat.at(0*2+1), "xgemm, x=d");
    formic::test::equal(8.7, mat.at(1*2+1), "xgemm, x=d");
  }

  // zgemm
  {
    const int n = 2;
    std::vector<std::complex<double> > mat(n*n, std::complex<double>(1.0, 0.4));
    const std::complex<double> alpha(0.5, 0.6);
    const std::complex<double> beta(0.2, -1.3);
    formic::zgemm('N', 'N', n, n, n, alpha, &mc0.at(0), n, &mc1.at(0), n, beta, &mat.at(0), n);
    formic::test::equal(std::complex<double>( 3.22,  1.78), mat.at(0*2+0), "zgemm");
    formic::test::equal(std::complex<double>(-4.68,  3.28), mat.at(1*2+0), "zgemm");
    formic::test::equal(std::complex<double>( 4.32, -4.22), mat.at(0*2+1), "zgemm");
    formic::test::equal(std::complex<double>( 4.72,  3.58), mat.at(1*2+1), "zgemm");
  }

  // xgemm, x=z
  {
    const int n = 2;
    std::vector<std::complex<double> > mat(n*n, std::complex<double>(1.0, 0.4));
    const std::complex<double> alpha(0.5, 0.6);
    const std::complex<double> beta(0.2, -1.3);
    formic::xgemm('N', 'N', n, n, n, alpha, &mc0.at(0), n, &mc1.at(0), n, beta, &mat.at(0), n);
    formic::test::equal(std::complex<double>( 3.22,  1.78), mat.at(0*2+0), "xgemm, x=z");
    formic::test::equal(std::complex<double>(-4.68,  3.28), mat.at(1*2+0), "xgemm, x=z");
    formic::test::equal(std::complex<double>( 4.32, -4.22), mat.at(0*2+1), "xgemm, x=z");
    formic::test::equal(std::complex<double>( 4.72,  3.58), mat.at(1*2+1), "xgemm, x=z");
  }

  // dsyev
  {
    const int n = 2;
    const int lwork = 3*n;
    std::vector<double> mat(md1);
    std::vector<double> w(n, 0.0);
    std::vector<double> work(lwork, 0.0);
    int info = 1;
    formic::dsyev('V', 'U', n, &mat.at(0), n, &w.at(0), &work.at(0), lwork, info);
    if (info != 0)
      throw std::runtime_error((boost::format("dsyev failed with info = %i") % info).str());
    formic::test::equal(w.at(0), 2.0 - std::sqrt(5.0), "dsyev");
    formic::test::equal(w.at(1), 2.0 + std::sqrt(5.0), "dsyev");
    formic::test::equal(md1.at(0*2+0)*mat.at(0*2+0) + md1.at(1*2+0)*mat.at(0*2+1), w.at(0)*mat.at(0*2+0), "dsyev");
    formic::test::equal(md1.at(0*2+1)*mat.at(0*2+0) + md1.at(1*2+1)*mat.at(0*2+1), w.at(0)*mat.at(0*2+1), "dsyev");
    formic::test::equal(md1.at(0*2+0)*mat.at(1*2+0) + md1.at(1*2+0)*mat.at(1*2+1), w.at(1)*mat.at(1*2+0), "dsyev");
    formic::test::equal(md1.at(0*2+1)*mat.at(1*2+0) + md1.at(1*2+1)*mat.at(1*2+1), w.at(1)*mat.at(1*2+1), "dsyev");
  }

  // xsyev, x=d
  {
    const int n = 2;
    const int lwork = 3*n;
    std::vector<double> mat(md1);
    std::vector<double> w(n, 0.0);
    std::vector<double> work(lwork, 0.0);
    std::vector<double> rwork(3*n, 0.0);
    int info = 1;
    formic::xsyev('V', 'U', n, &mat.at(0), n, &w.at(0), &work.at(0), lwork, &rwork.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("xsyev, x=d failed with info = %i") % info).str());
    formic::test::equal(w.at(0), 2.0 - std::sqrt(5.0), "xsyev, x=d");
    formic::test::equal(w.at(1), 2.0 + std::sqrt(5.0), "xsyev, x=d");
    formic::test::equal(md1.at(0*2+0)*mat.at(0*2+0) + md1.at(1*2+0)*mat.at(0*2+1), w.at(0)*mat.at(0*2+0), "xsyev, x=d");
    formic::test::equal(md1.at(0*2+1)*mat.at(0*2+0) + md1.at(1*2+1)*mat.at(0*2+1), w.at(0)*mat.at(0*2+1), "xsyev, x=d");
    formic::test::equal(md1.at(0*2+0)*mat.at(1*2+0) + md1.at(1*2+0)*mat.at(1*2+1), w.at(1)*mat.at(1*2+0), "xsyev, x=d");
    formic::test::equal(md1.at(0*2+1)*mat.at(1*2+0) + md1.at(1*2+1)*mat.at(1*2+1), w.at(1)*mat.at(1*2+1), "xsyev, x=d");
  }

  // zheev
  {
    const int n = 2;
    const int lwork = 2*n;
    std::vector<std::complex<double> > mat(mc1);
    std::vector<double> w(n, 0.0);
    std::vector<std::complex<double> > work(lwork, std::complex<double>(0.0,0.0));
    std::vector<double> rwork(3*n, 0.0);
    int info = 1;
    formic::zheev('V', 'U', n, &mat.at(0), n, &w.at(0), &work.at(0), lwork, &rwork.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("zheev failed with info = %i") % info).str());
    formic::test::equal(w.at(0), 2.0 - std::sqrt(5.0), "zheev");
    formic::test::equal(w.at(1), 2.0 + std::sqrt(5.0), "zheev");
    formic::test::equal(mc1.at(0*2+0)*mat.at(0*2+0) + mc1.at(1*2+0)*mat.at(0*2+1), w.at(0)*mat.at(0*2+0), "zheev");
    formic::test::equal(mc1.at(0*2+1)*mat.at(0*2+0) + mc1.at(1*2+1)*mat.at(0*2+1), w.at(0)*mat.at(0*2+1), "zheev");
    formic::test::equal(mc1.at(0*2+0)*mat.at(1*2+0) + mc1.at(1*2+0)*mat.at(1*2+1), w.at(1)*mat.at(1*2+0), "zheev");
    formic::test::equal(mc1.at(0*2+1)*mat.at(1*2+0) + mc1.at(1*2+1)*mat.at(1*2+1), w.at(1)*mat.at(1*2+1), "zheev");
  }

  // xsyev, x=z
  {
    const int n = 2;
    const int lwork = 2*n;
    std::vector<std::complex<double> > mat(mc1);
    std::vector<double> w(n, 0.0);
    std::vector<std::complex<double> > work(lwork, std::complex<double>(0.0,0.0));
    std::vector<double> rwork(3*n, 0.0);
    int info = 1;
    formic::xsyev('V', 'U', n, &mat.at(0), n, &w.at(0), &work.at(0), lwork, &rwork.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("xsyev, x=z failed with info = %i") % info).str());
    formic::test::equal(w.at(0), 2.0 - std::sqrt(5.0), "xsyev, x=z");
    formic::test::equal(w.at(1), 2.0 + std::sqrt(5.0), "xsyev, x=z");
    formic::test::equal(mc1.at(0*2+0)*mat.at(0*2+0) + mc1.at(1*2+0)*mat.at(0*2+1), w.at(0)*mat.at(0*2+0), "xsyev, x=z");
    formic::test::equal(mc1.at(0*2+1)*mat.at(0*2+0) + mc1.at(1*2+1)*mat.at(0*2+1), w.at(0)*mat.at(0*2+1), "xsyev, x=z");
    formic::test::equal(mc1.at(0*2+0)*mat.at(1*2+0) + mc1.at(1*2+0)*mat.at(1*2+1), w.at(1)*mat.at(1*2+0), "xsyev, x=z");
    formic::test::equal(mc1.at(0*2+1)*mat.at(1*2+0) + mc1.at(1*2+1)*mat.at(1*2+1), w.at(1)*mat.at(1*2+1), "xsyev, x=z");
  }

  // dsygv
  {
    const int n = 2;
    const int lwork = 3*n;
    std::vector<double> mata(md1);
    std::vector<double> matb(md2);
    std::vector<double> w(n, 0.0);
    std::vector<double> work(lwork, 0.0);
    int info = 1;
    formic::dsygv(1, 'V', 'U', n, &mata.at(0), n, &matb.at(0), n, &w.at(0), &work.at(0), lwork, info);
    if (info != 0)
      throw std::runtime_error((boost::format("dsygv failed with info = %i") % info).str());
    formic::test::equal(         md1.at(0*2+0)*mata.at(0*2+0) + md1.at(1*2+0)*mata.at(0*2+1),
                        w.at(0)*(md2.at(0*2+0)*mata.at(0*2+0) + md2.at(1*2+0)*mata.at(0*2+1)), "dsygv");
    formic::test::equal(         md1.at(0*2+1)*mata.at(0*2+0) + md1.at(1*2+1)*mata.at(0*2+1),
                        w.at(0)*(md2.at(0*2+1)*mata.at(0*2+0) + md2.at(1*2+1)*mata.at(0*2+1)), "dsygv");
    formic::test::equal(         md1.at(0*2+0)*mata.at(1*2+0) + md1.at(1*2+0)*mata.at(1*2+1),
                        w.at(1)*(md2.at(0*2+0)*mata.at(1*2+0) + md2.at(1*2+0)*mata.at(1*2+1)), "dsygv");
    formic::test::equal(         md1.at(0*2+1)*mata.at(1*2+0) + md1.at(1*2+1)*mata.at(1*2+1),
                        w.at(1)*(md2.at(0*2+1)*mata.at(1*2+0) + md2.at(1*2+1)*mata.at(1*2+1)), "dsygv");
  }

  // xsygv, x=d
  {
    const int n = 2;
    const int lwork = 3*n;
    std::vector<double> mata(md1);
    std::vector<double> matb(md2);
    std::vector<double> w(n, 0.0);
    std::vector<double> work(lwork, 0.0);
    std::vector<double> rwork(3*n, 0.0);
    int info = 1;
    formic::xsygv(1, 'V', 'U', n, &mata.at(0), n, &matb.at(0), n, &w.at(0), &work.at(0), lwork, &rwork.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("xsygv, x=d failed with info = %i") % info).str());
    formic::test::equal(         md1.at(0*2+0)*mata.at(0*2+0) + md1.at(1*2+0)*mata.at(0*2+1),
                        w.at(0)*(md2.at(0*2+0)*mata.at(0*2+0) + md2.at(1*2+0)*mata.at(0*2+1)), "xsygv, x=d");
    formic::test::equal(         md1.at(0*2+1)*mata.at(0*2+0) + md1.at(1*2+1)*mata.at(0*2+1),
                        w.at(0)*(md2.at(0*2+1)*mata.at(0*2+0) + md2.at(1*2+1)*mata.at(0*2+1)), "xsygv, x=d");
    formic::test::equal(         md1.at(0*2+0)*mata.at(1*2+0) + md1.at(1*2+0)*mata.at(1*2+1),
                        w.at(1)*(md2.at(0*2+0)*mata.at(1*2+0) + md2.at(1*2+0)*mata.at(1*2+1)), "xsygv, x=d");
    formic::test::equal(         md1.at(0*2+1)*mata.at(1*2+0) + md1.at(1*2+1)*mata.at(1*2+1),
                        w.at(1)*(md2.at(0*2+1)*mata.at(1*2+0) + md2.at(1*2+1)*mata.at(1*2+1)), "xsygv, x=d");
  }

  // zhegv
  {
    const int n = 2;
    const int lwork = 2*n;
    std::vector<std::complex<double> > mata(mc1);
    std::vector<std::complex<double> > matb(mc2);
    std::vector<double> w(n, 0.0);
    std::vector<std::complex<double> > work(lwork, std::complex<double>(0.0,0.0));
    std::vector<double> rwork(3*n, 0.0);
    int info = 1;
    formic::zhegv(1, 'V', 'U', n, &mata.at(0), n, &matb.at(0), n, &w.at(0), &work.at(0), lwork, &rwork.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("zhegv failed with info = %i") % info).str());
    formic::test::equal(         mc1.at(0*2+0)*mata.at(0*2+0) + mc1.at(1*2+0)*mata.at(0*2+1),
                        w.at(0)*(mc2.at(0*2+0)*mata.at(0*2+0) + mc2.at(1*2+0)*mata.at(0*2+1)), "zhegv");
    formic::test::equal(         mc1.at(0*2+1)*mata.at(0*2+0) + mc1.at(1*2+1)*mata.at(0*2+1),
                        w.at(0)*(mc2.at(0*2+1)*mata.at(0*2+0) + mc2.at(1*2+1)*mata.at(0*2+1)), "zhegv");
    formic::test::equal(         mc1.at(0*2+0)*mata.at(1*2+0) + mc1.at(1*2+0)*mata.at(1*2+1),
                        w.at(1)*(mc2.at(0*2+0)*mata.at(1*2+0) + mc2.at(1*2+0)*mata.at(1*2+1)), "zhegv");
    formic::test::equal(         mc1.at(0*2+1)*mata.at(1*2+0) + mc1.at(1*2+1)*mata.at(1*2+1),
                        w.at(1)*(mc2.at(0*2+1)*mata.at(1*2+0) + mc2.at(1*2+1)*mata.at(1*2+1)), "zhegv");
  }

  // xsygv, x=z
  {
    const int n = 2;
    const int lwork = 2*n;
    std::vector<std::complex<double> > mata(mc1);
    std::vector<std::complex<double> > matb(mc2);
    std::vector<double> w(n, 0.0);
    std::vector<std::complex<double> > work(lwork, std::complex<double>(0.0,0.0));
    std::vector<double> rwork(3*n, 0.0);
    int info = 1;
    formic::xsygv(1, 'V', 'U', n, &mata.at(0), n, &matb.at(0), n, &w.at(0), &work.at(0), lwork, &rwork.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("xsygv, x=z failed with info = %i") % info).str());
    formic::test::equal(         mc1.at(0*2+0)*mata.at(0*2+0) + mc1.at(1*2+0)*mata.at(0*2+1),
                        w.at(0)*(mc2.at(0*2+0)*mata.at(0*2+0) + mc2.at(1*2+0)*mata.at(0*2+1)), "xsygv, x=z");
    formic::test::equal(         mc1.at(0*2+1)*mata.at(0*2+0) + mc1.at(1*2+1)*mata.at(0*2+1),
                        w.at(0)*(mc2.at(0*2+1)*mata.at(0*2+0) + mc2.at(1*2+1)*mata.at(0*2+1)), "xsygv, x=z");
    formic::test::equal(         mc1.at(0*2+0)*mata.at(1*2+0) + mc1.at(1*2+0)*mata.at(1*2+1),
                        w.at(1)*(mc2.at(0*2+0)*mata.at(1*2+0) + mc2.at(1*2+0)*mata.at(1*2+1)), "xsygv, x=z");
    formic::test::equal(         mc1.at(0*2+1)*mata.at(1*2+0) + mc1.at(1*2+1)*mata.at(1*2+1),
                        w.at(1)*(mc2.at(0*2+1)*mata.at(1*2+0) + mc2.at(1*2+1)*mata.at(1*2+1)), "xsygv, x=z");
  }

  // dgeev
  {
    const int n = 2;
    const int lwork = 4*n;
    std::vector<double> mat(md0);
    std::vector<double> vl(md0);
    std::vector<double> vr(md0);
    std::vector<double> wr(n, 7.0);
    std::vector<double> wi(n, 7.0);
    std::vector<double> work(lwork, 0.0);
    int info = 1;
    formic::dgeev('V', 'V', n, &mat.at(0), n, &wr.at(0), &wi.at(0), &vl.at(0), n, &vr.at(0), n, &work.at(0), lwork, info);
    if (info != 0)
      throw std::runtime_error((boost::format("dgeev failed with info = %i") % info).str());
    formic::test::equal(wi.at(0), 0.0, "dgeev");
    formic::test::equal(wi.at(1), 0.0, "dgeev");
    formic::test::equal(md0.at(0*2+0)*vr.at(0*2+0) + md0.at(1*2+0)*vr.at(0*2+1), wr.at(0)*vr.at(0*2+0), "dgeev");
    formic::test::equal(md0.at(0*2+1)*vr.at(0*2+0) + md0.at(1*2+1)*vr.at(0*2+1), wr.at(0)*vr.at(0*2+1), "dgeev");
    formic::test::equal(md0.at(0*2+0)*vr.at(1*2+0) + md0.at(1*2+0)*vr.at(1*2+1), wr.at(1)*vr.at(1*2+0), "dgeev");
    formic::test::equal(md0.at(0*2+1)*vr.at(1*2+0) + md0.at(1*2+1)*vr.at(1*2+1), wr.at(1)*vr.at(1*2+1), "dgeev");
  }

  // xgeev, x=d
  {
    const int n = 2;
    const int lwork = 4*n;
    std::vector<double> mat(md0);
    std::vector<double> vl(md0);
    std::vector<double> vr(md0);
    std::vector<double> w(n, 7.0);
    std::vector<double> wr(n, 7.0);
    std::vector<double> wi(n, 7.0);
    std::vector<double> work(lwork, 0.0);
    std::vector<double> rwork(2*n, 0.0);
    int info = 1;
    formic::xgeev('V', 'V', n, &mat.at(0), n, &w.at(0), &wr.at(0), &wi.at(0), &vl.at(0), n, &vr.at(0), n,
                  &work.at(0), lwork, &rwork.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("xgeev, x=d failed with info = %i") % info).str());
    formic::test::equal(wi.at(0), 0.0, "xgeev, x=d");
    formic::test::equal(wi.at(1), 0.0, "xgeev, x=d");
    formic::test::equal(md0.at(0*2+0)*vr.at(0*2+0) + md0.at(1*2+0)*vr.at(0*2+1), wr.at(0)*vr.at(0*2+0), "xgeev, x=d");
    formic::test::equal(md0.at(0*2+1)*vr.at(0*2+0) + md0.at(1*2+1)*vr.at(0*2+1), wr.at(0)*vr.at(0*2+1), "xgeev, x=d");
    formic::test::equal(md0.at(0*2+0)*vr.at(1*2+0) + md0.at(1*2+0)*vr.at(1*2+1), wr.at(1)*vr.at(1*2+0), "xgeev, x=d");
    formic::test::equal(md0.at(0*2+1)*vr.at(1*2+0) + md0.at(1*2+1)*vr.at(1*2+1), wr.at(1)*vr.at(1*2+1), "xgeev, x=d");
  }

  // zgeev
  {
    const int n = 2;
    const int lwork = 2*n;
    std::vector<std::complex<double> > mat(mc0);
    std::vector<std::complex<double> > vl(mc0);
    std::vector<std::complex<double> > vr(mc0);
    std::vector<std::complex<double> > w(n, std::complex<double>(7.0,1.0));
    std::vector<std::complex<double> > work(lwork, std::complex<double>(7.0,1.0));
    std::vector<double> rwork(2*n, 0.0);
    int info = 1;
    formic::zgeev('V', 'V', n, &mat.at(0), n, &w.at(0), &vl.at(0), n, &vr.at(0), n, &work.at(0), lwork, &rwork.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("zgeev failed with info = %i") % info).str());
    const std::complex<double> e0(3.5,  std::sqrt(7.0/4.0));
    const std::complex<double> e1(3.5, -std::sqrt(7.0/4.0));
    if (w.at(0).imag() > 0.0) {
      formic::test::equal(w.at(0), e0, "zgeev");
      formic::test::equal(w.at(1), e1, "zgeev");
    } else {
      formic::test::equal(w.at(0), e1, "zgeev");
      formic::test::equal(w.at(1), e0, "zgeev");
    }
    formic::test::equal(mc0.at(0*2+0)*vr.at(0*2+0) + mc0.at(1*2+0)*vr.at(0*2+1), w.at(0)*vr.at(0*2+0), "zgeev");
    formic::test::equal(mc0.at(0*2+1)*vr.at(0*2+0) + mc0.at(1*2+1)*vr.at(0*2+1), w.at(0)*vr.at(0*2+1), "zgeev");
    formic::test::equal(mc0.at(0*2+0)*vr.at(1*2+0) + mc0.at(1*2+0)*vr.at(1*2+1), w.at(1)*vr.at(1*2+0), "zgeev");
    formic::test::equal(mc0.at(0*2+1)*vr.at(1*2+0) + mc0.at(1*2+1)*vr.at(1*2+1), w.at(1)*vr.at(1*2+1), "zgeev");
  }

  // xgeev, x=z
  {
    const int n = 2;
    const int lwork = 2*n;
    std::vector<std::complex<double> > mat(mc0);
    std::vector<std::complex<double> > vl(mc0);
    std::vector<std::complex<double> > vr(mc0);
    std::vector<std::complex<double> > w(n, std::complex<double>(7.0,1.0));
    std::vector<std::complex<double> > work(lwork, std::complex<double>(7.0,1.0));
    std::vector<double> rwork(2*n, 0.0);
    std::vector<double> wr(n, 4.0);
    std::vector<double> wi(n, 4.0);
    int info = 1;
    formic::xgeev('V', 'V', n, &mat.at(0), n, &w.at(0), &wr.at(0), &wi.at(0), &vl.at(0), n, &vr.at(0), n,
                  &work.at(0), lwork, &rwork.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("xgeev, x=z failed with info = %i") % info).str());
    const std::complex<double> e0(3.5,  std::sqrt(7.0/4.0));
    const std::complex<double> e1(3.5, -std::sqrt(7.0/4.0));
    formic::test::equal(wr.at(0), 3.5, "xgeev, x=z");
    formic::test::equal(wr.at(1), 3.5, "xgeev, x=z");
    if (w.at(0).imag() > 0.0) {
      formic::test::equal(w.at(0), e0, "xgeev, x=z");
      formic::test::equal(w.at(1), e1, "xgeev, x=z");
      formic::test::equal(wi.at(0),  std::sqrt(7.0/4.0), "xgeev, x=z");
      formic::test::equal(wi.at(1), -std::sqrt(7.0/4.0), "xgeev, x=z");
    } else {
      formic::test::equal(w.at(0), e1, "xgeev, x=z");
      formic::test::equal(w.at(1), e0, "xgeev, x=z");
      formic::test::equal(wi.at(0), -std::sqrt(7.0/4.0), "xgeev, x=z");
      formic::test::equal(wi.at(1),  std::sqrt(7.0/4.0), "xgeev, x=z");
    }
    formic::test::equal(mc0.at(0*2+0)*vr.at(0*2+0) + mc0.at(1*2+0)*vr.at(0*2+1), w.at(0)*vr.at(0*2+0), "xgeev, x=z");
    formic::test::equal(mc0.at(0*2+1)*vr.at(0*2+0) + mc0.at(1*2+1)*vr.at(0*2+1), w.at(0)*vr.at(0*2+1), "xgeev, x=z");
    formic::test::equal(mc0.at(0*2+0)*vr.at(1*2+0) + mc0.at(1*2+0)*vr.at(1*2+1), w.at(1)*vr.at(1*2+0), "xgeev, x=z");
    formic::test::equal(mc0.at(0*2+1)*vr.at(1*2+0) + mc0.at(1*2+1)*vr.at(1*2+1), w.at(1)*vr.at(1*2+1), "xgeev, x=z");
  }

  // dgesvd
  {
    const int n = 2;
    const int lwork = 5*n;
    std::vector<double> mat(n*n, 0.0);
    mat.at(0*2+0) = 3.0;
    mat.at(1*2+0) = 6.0;
    mat.at(0*2+1) = 5.0;
    mat.at(1*2+1) = 9.0;
    std::vector<double> mata(mat);
    std::vector<double> matu(n*n, 0.0);
    std::vector<double> matvt(n*n, 0.0);
    std::vector<double> sigma(n, 0.0);
    std::vector<double> work(lwork, 0.0);
    int info = 1;
    formic::dgesvd('A', 'A', n, n, &mata.at(0), n, &sigma.at(0), &matu.at(0), n, &matvt.at(0), n, &work.at(0), lwork, info);
    if (info != 0)
      throw std::runtime_error((boost::format("dgesvd failed with info = %i") % info).str());
    matvt.at(0*2+0) *= sigma.at(0);
    matvt.at(1*2+0) *= sigma.at(0);
    matvt.at(0*2+1) *= sigma.at(1);
    matvt.at(1*2+1) *= sigma.at(1);
    formic::dgemm('N', 'N', n, n, n, 1.0, &matu.at(0), n, &matvt.at(0), n, 0.0, &mata.at(0), n);
    formic::test::equal(mat, mata, "dgesvd");
  }

  // xgesvd, x=d
  {
    const int n = 2;
    const int lwork = 5*n;
    std::vector<double> mat(n*n, 0.0);
    mat.at(0*2+0) = 3.0;
    mat.at(1*2+0) = 6.0;
    mat.at(0*2+1) = 5.0;
    mat.at(1*2+1) = 9.0;
    std::vector<double> mata(mat);
    std::vector<double> matu(n*n, 0.0);
    std::vector<double> matvt(n*n, 0.0);
    std::vector<double> sigma(n, 0.0);
    std::vector<double> work(lwork, 0.0);
    std::vector<double> rwork(5*n, 0.0);
    int info = 1;
    formic::xgesvd('A', 'A', n, n, &mata.at(0), n, &sigma.at(0), &matu.at(0), n, &matvt.at(0), n, &work.at(0), lwork, &rwork.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("xgesvd, x=d failed with info = %i") % info).str());
    matvt.at(0*2+0) *= sigma.at(0);
    matvt.at(1*2+0) *= sigma.at(0);
    matvt.at(0*2+1) *= sigma.at(1);
    matvt.at(1*2+1) *= sigma.at(1);
    formic::xgemm('N', 'N', n, n, n, 1.0, &matu.at(0), n, &matvt.at(0), n, 0.0, &mata.at(0), n);
    formic::test::equal(mat, mata, "xgesvd, x=d");
  }

  // zgesvd
  {
    const int n = 2;
    const int lwork = 3*n;
    std::vector<std::complex<double> > mat(n*n, 0.0);
    mat.at(0*2+0) = std::complex<double>(3.0, 0.2);
    mat.at(1*2+0) = std::complex<double>(6.0, 0.1);
    mat.at(0*2+1) = std::complex<double>(5.0,-0.1);
    mat.at(1*2+1) = std::complex<double>(9.0,-0.2);
    std::vector<std::complex<double> > mata(mat);
    std::vector<std::complex<double> > matu(n*n, 0.0);
    std::vector<std::complex<double> > matvt(n*n, 0.0);
    std::vector<double> sigma(n, 0.0);
    std::vector<std::complex<double> > work(lwork, 0.0);
    std::vector<double> rwork(5*n, 0.0);
    int info = 1;
    formic::zgesvd('A', 'A', n, n, &mata.at(0), n, &sigma.at(0), &matu.at(0), n, &matvt.at(0), n, &work.at(0), lwork, &rwork.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("zgesvd failed with info = %i") % info).str());
    matvt.at(0*2+0) *= sigma.at(0);
    matvt.at(1*2+0) *= sigma.at(0);
    matvt.at(0*2+1) *= sigma.at(1);
    matvt.at(1*2+1) *= sigma.at(1);
    formic::zgemm('N', 'N', n, n, n, std::complex<double>(1.0, 0.0), &matu.at(0), n, &matvt.at(0), n,
                  std::complex<double>(0.0, 0.0), &mata.at(0), n);
    formic::test::equal(mat, mata, "zgesvd");
  }

  // xgesvd, x=z
  {
    const int n = 2;
    const int lwork = 3*n;
    std::vector<std::complex<double> > mat(n*n, 0.0);
    mat.at(0*2+0) = std::complex<double>(3.0, 0.2);
    mat.at(1*2+0) = std::complex<double>(6.0, 0.1);
    mat.at(0*2+1) = std::complex<double>(5.0,-0.1);
    mat.at(1*2+1) = std::complex<double>(9.0,-0.2);
    std::vector<std::complex<double> > mata(mat);
    std::vector<std::complex<double> > matu(n*n, 0.0);
    std::vector<std::complex<double> > matvt(n*n, 0.0);
    std::vector<double> sigma(n, 0.0);
    std::vector<std::complex<double> > work(lwork, 0.0);
    std::vector<double> rwork(5*n, 0.0);
    int info = 1;
    formic::xgesvd('A', 'A', n, n, &mata.at(0), n, &sigma.at(0), &matu.at(0), n, &matvt.at(0), n, &work.at(0), lwork, &rwork.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("xgesvd, x=z failed with info = %i") % info).str());
    matvt.at(0*2+0) *= sigma.at(0);
    matvt.at(1*2+0) *= sigma.at(0);
    matvt.at(0*2+1) *= sigma.at(1);
    matvt.at(1*2+1) *= sigma.at(1);
    formic::xgemm('N', 'N', n, n, n, std::complex<double>(1.0, 0.0), &matu.at(0), n, &matvt.at(0), n,
                  std::complex<double>(0.0, 0.0), &mata.at(0), n);
    formic::test::equal(mat, mata, "xgesvd, x=z");
  }

  // dgetrf
  {
    const int n = 2;
    std::vector<double> mat(n*n, 0.0);
    mat.at(0*2+0) = 3.0; //3.0;
    mat.at(1*2+0) = 6.0; //6.0;
    mat.at(0*2+1) = 5.0;
    mat.at(1*2+1) = 9.0; //9.0;
    std::vector<double> matl(mat);
    std::vector<double> matu(mat);
    std::vector<double> mata(mat);
    std::vector<int> ipiv(n, 0);
    int info = 1;
    formic::dgetrf(n, n, &matu.at(0), n, &ipiv.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("dgetrf failed with info = %i") % info).str());
    matl.at(0*2+0) = 1.0;
    matl.at(1*2+1) = 1.0;
    matl.at(1*2+0) = 0.0;
    matl.at(0*2+1) = matu.at(0*2+1);
    matu.at(0*2+1) = 0.0;
    formic::dgemm('N', 'N', n, n, n, 1.0, &matl.at(0), n, &matu.at(0), n, 0.0, &mata.at(0), n);
    if (ipiv.at(0) == 2) {
      std::swap(mata.at(0*2+0), mata.at(0*2+1));
      std::swap(mata.at(1*2+0), mata.at(1*2+1));
    }
    if (ipiv.at(1) == 1) {
      std::swap(mata.at(0*2+0), mata.at(0*2+1));
      std::swap(mata.at(1*2+0), mata.at(1*2+1));
    }
    formic::test::equal(mat, mata, "dgetrf");
  }

  // xgetrf, x=d
  {
    const int n = 2;
    std::vector<double> mat(n*n, 0.0);
    mat.at(0*2+0) = 3.0; //3.0;
    mat.at(1*2+0) = 6.0; //6.0;
    mat.at(0*2+1) = 5.0;
    mat.at(1*2+1) = 9.0; //9.0;
    std::vector<double> matl(mat);
    std::vector<double> matu(mat);
    std::vector<double> mata(mat);
    std::vector<int> ipiv(n, 0);
    int info = 1;
    formic::xgetrf(n, n, &matu.at(0), n, &ipiv.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("xgetrf, x=d failed with info = %i") % info).str());
    matl.at(0*2+0) = 1.0;
    matl.at(1*2+1) = 1.0;
    matl.at(1*2+0) = 0.0;
    matl.at(0*2+1) = matu.at(0*2+1);
    matu.at(0*2+1) = 0.0;
    formic::dgemm('N', 'N', n, n, n, 1.0, &matl.at(0), n, &matu.at(0), n, 0.0, &mata.at(0), n);
    if (ipiv.at(0) == 2) {
      std::swap(mata.at(0*2+0), mata.at(0*2+1));
      std::swap(mata.at(1*2+0), mata.at(1*2+1));
    }
    if (ipiv.at(1) == 1) {
      std::swap(mata.at(0*2+0), mata.at(0*2+1));
      std::swap(mata.at(1*2+0), mata.at(1*2+1));
    }
    formic::test::equal(mat, mata, "xgetrf, x=d");
  }

  // zgetrf
  {
    const int n = 2;
    std::vector<std::complex<double> > mat(n*n);
    mat.at(0*2+0) = std::complex<double>(3.0,  0.1);
    mat.at(1*2+0) = std::complex<double>(6.0,  0.2);
    mat.at(0*2+1) = std::complex<double>(5.0, -0.1);
    mat.at(1*2+1) = std::complex<double>(9.0, -0.2);
    std::vector<std::complex<double> > matl(mat);
    std::vector<std::complex<double> > matu(mat);
    std::vector<std::complex<double> > mata(mat);
    std::vector<int> ipiv(n, 0);
    int info = 1;
    formic::zgetrf(n, n, &matu.at(0), n, &ipiv.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("zgetrf failed with info = %i") % info).str());
    matl.at(0*2+0) = std::complex<double>(1.0, 0.0);
    matl.at(1*2+1) = std::complex<double>(1.0, 0.0);
    matl.at(1*2+0) = std::complex<double>(0.0, 0.0);
    matl.at(0*2+1) = matu.at(0*2+1);
    matu.at(0*2+1) = std::complex<double>(0.0, 0.0);
    formic::zgemm('N', 'N', n, n, n, std::complex<double>(1.0,0.0), &matl.at(0), n, &matu.at(0), n, std::complex<double>(0.0,0.0), &mata.at(0), n);
    if (ipiv.at(0) == 2) {
      std::swap(mata.at(0*2+0), mata.at(0*2+1));
      std::swap(mata.at(1*2+0), mata.at(1*2+1));
    }
    if (ipiv.at(1) == 1) {
      std::swap(mata.at(0*2+0), mata.at(0*2+1));
      std::swap(mata.at(1*2+0), mata.at(1*2+1));
    }
    formic::test::equal(mat, mata, "zgetrf");
  }

  // xgetrf, x=z
  {
    const int n = 2;
    std::vector<std::complex<double> > mat(n*n);
    mat.at(0*2+0) = std::complex<double>(3.0,  0.1);
    mat.at(1*2+0) = std::complex<double>(6.0,  0.2);
    mat.at(0*2+1) = std::complex<double>(5.0, -0.1);
    mat.at(1*2+1) = std::complex<double>(9.0, -0.2);
    std::vector<std::complex<double> > matl(mat);
    std::vector<std::complex<double> > matu(mat);
    std::vector<std::complex<double> > mata(mat);
    std::vector<int> ipiv(n, 0);
    int info = 1;
    formic::xgetrf(n, n, &matu.at(0), n, &ipiv.at(0), info);
    if (info != 0)
      throw std::runtime_error((boost::format("xgetrf, x=z failed with info = %i") % info).str());
    matl.at(0*2+0) = std::complex<double>(1.0, 0.0);
    matl.at(1*2+1) = std::complex<double>(1.0, 0.0);
    matl.at(1*2+0) = std::complex<double>(0.0, 0.0);
    matl.at(0*2+1) = matu.at(0*2+1);
    matu.at(0*2+1) = std::complex<double>(0.0, 0.0);
    formic::zgemm('N', 'N', n, n, n, std::complex<double>(1.0,0.0), &matl.at(0), n, &matu.at(0), n, std::complex<double>(0.0,0.0), &mata.at(0), n);
    if (ipiv.at(0) == 2) {
      std::swap(mata.at(0*2+0), mata.at(0*2+1));
      std::swap(mata.at(1*2+0), mata.at(1*2+1));
    }
    if (ipiv.at(1) == 1) {
      std::swap(mata.at(0*2+0), mata.at(0*2+1));
      std::swap(mata.at(1*2+0), mata.at(1*2+1));
    }
    formic::test::equal(mat, mata, "xgetrf, x=z");
  }

}

int main() {

  try {

    // set the ouput file to "not finished"
    {
      std::ofstream ofs("lapack_test_output.txt", std::ios_base::out | std::ios_base::trunc);
      if (!ofs.is_open())
        throw std::runtime_error("failed to open output file");
      ofs << int(1) << " " << std::endl;
      ofs.close();
    }

    // run the tests
    ::formic_lapack_tests();

    // set the output file to "finished"
    {
      std::ofstream ofs("lapack_test_output.txt", std::ios_base::out | std::ios_base::trunc);
      if (!ofs.is_open())
        throw std::runtime_error("failed to open output file");
      ofs << int(2) << " " << std::endl;
      ofs.close();
    }

  } catch (std::exception &e) {
    std::cout << "Exception occurred:  " << e.what() << std::endl;
    return EXIT_FAILURE;
  } catch (...) {
    std::cout << "Exception occurred." << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;

}
